%% 灰色预测模型汇总
clear,clc
% 生成各类型序列
x=0:0.1:10;
f1=diff(exp(-0.02*x)+1);% % 指数下降类序列，适合GM(1,1)预测
f2=diff(1e2*exp(x)-1e1*exp(0.001*x)+2);% 双指数综合，可选用GM(2,1)预测
f3=diff(2.7*exp(-0.4*x)+4*x+1);% 指数与一次综合，可选用DGM(2,1)预测
f4=diff(1./(1+50*exp(-0.01*x)));% s型曲线，可选用Verhulst预测
f=[f1;f2;f3;f4];
% 预测数据，与原始数据作比较
for i=1:4
    data=f(i,:);
    n=length(data);% 原始数据长度
    data0=diff(data);
    data1=cumsum(data);
    z_data=1/2*(data1(2:end)+data1(1:end-1));
    if(i==1)% GM(1,1)
        B=[-z_data(1:end)',ones(n-1,1)];
        Y=[data(2:end)]';
        u=(B'*B)\B'*Y;
        pre=@(k)(data(1)-u(2)/u(1))*exp(-u(1)*(k-1))+u(2)/u(1);
        fp=pre(1:n);fp=diff(fp);
    end
    if(i==2)% GM(2,1)
        B=[-data(2:end)',-z_data(1:end)',ones(n-1,1)];
        Y=data0';
        u=(B'*B)\B'*Y;
        syms x_un(t)
        x_un=dsolve(diff(x_un,2)+u(1)*diff(x_un)+u(2)*x_un==u(3), ...
            x_un(0)==data1(1),x_un(n-1)==data1(end));% 求解白化方程
        xt=vpa(x_un,6);
        fp=double(subs(x_un,t,0:n-1));
    end
    if(i==3)% DGM(2,1)
        B=[-data(2:end)',ones(n-1,1)];
        Y=data0';
        u=(B'*B)\B'*Y;
        pre=@(k)(u(2)/u(1)^2-data(1)/u(1))*exp(-u(1)*(k-1))+...
        u(2)/u(1)*(k-1)+(1+u(1))/u(1)*data(1)-u(2)/u(1)^2;
        fp=pre(1:n);fp=diff(fp);
    end
    if(i==4)% Verhulst
        B=[-z_data',z_data'.^2];
        Y=data(2:end)';
        u=(B'*B)\B'*Y;
        pre=@(k)u(1)*data(1)./(u(2)*data(1)+(u(1)-u(2)*data(1))*exp(u(1)*(k-1)));
        fp=pre(1:n);fp=diff(fp);
    end
    figure(i)
    plot(data,'b*-')
    hold on
    plot(fp,'r*-')
    hold off
    legend('实际值','预测值')
end